setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/20201216")
file=read.csv("at.least.two.ASH.in.DC.BF.1.add.eQTL.GWAS.add.DEG.csv",header = T)
snp=unique(as.character(file$avsnp150))
link1="http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/dlpfc/"
link2="http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/hippo/"

###http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/hippo/rs4558409
###http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/dlpfc/rs4558409

links1=paste0(link1,snp)
links2=paste0(link2,snp)
links=c(links1,links2)
write.table(links,"./20201222/links.for.807ASH.LIBD.GWAS.txt",quote = F,row.names = F,sep="\t")
            
sel=list.files(path = "E:/迅雷下载/20201222.LIBD.data.807ASH/rs",pattern="rs",full.name=T)
i=2	#因为第一个文件是空的
tmp=read.csv(sel[i],header = T)
tmp=dplyr::select(tmp,snpRsNum,feature_id,Type)
rt=tmp

for(i in 3:length(sel)){
  tmp=read.csv(sel[i],header = T)
  if(dim(tmp)[1]>0){
    tmp=dplyr::select(tmp,snpRsNum,feature_id,Type)
    rt=rbind(rt,tmp) 
  }
}

rsid=unique(rt$snpRsNum)
i=1
rt1=rt[rt$snpRsNum==rsid[i],]
feature.id=paste(unique(as.character(rt1$feature_id)),collapse = ":")
types=paste(unique(as.character(rt1$Type)),collapse = ":")
rt1=data.frame(avsnp150=rt1[1,1],feature.id=feature.id,types=types)
rts=rt1
for(i in 2:length(rsid)){
  rt1=rt[rt$snpRsNum==rsid[i],]
  feature.id=paste(unique(as.character(rt1$feature_id)),collapse = ":")
  types=paste(unique(as.character(rt1$Type)),collapse = ":")
  rt1=data.frame(avsnp150=rt1[1,1],feature.id=feature.id,types=types)
  rts=rbind(rts,rt1)
}
names(rts)=c("avsnp150","LIBD.eQTL.feature.id","LIBD.eQTL.type")
file1=merge(file,rts,by="avsnp150",all.x = T)
write.csv(file1,"at.least.two.ASH.in.DC.BF.10.add.eQTL.GWAS.add.DEG.LIBDGWAS.csv")

file2=file1[!is.na(file1$Other)&!is.na(file1$LIBD.eQTL.type),]
length(unique(file2$unitID))
file2=file2[file2$gwas.no.na.num>0,]
length(unique(file2$unitID))
file2=file2[file2$DEG.validate.num>0,]
length(unique(file2$unitID))
file2=file2[!duplicated(file2$avsnp150),]
write.csv(file2,"Coordinates.site.csv",quote=F,row.names = F)

###给53K个ASH查询eQTL

setwd("E:/5hmc_file/2_5hmc_yjp_bam/ASM/")
file=read.csv("20201120/at.least.one.AShM.in.DC.add.BF.beta0.add.CCHC.csv",head=T)
snp=unique(as.character(file$avsnp150))
link1="http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/dlpfc/"
link2="http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/hippo/"

###http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/hippo/rs4558409
###http://eqtl.brainseq.org/db/ng-eqtl/snprsnum/csv/dlpfc/rs4558409

links1=paste0(link1,snp)
links2=paste0(link2,snp)
links=c(links1,links2)

write.table(links1,"./20210311/links.for.53K.ASH.LIBD.dlpfc.txt",quote = F,row.names = F,sep="\t")
write.table(links2,"./20210311/links.for.53K.ASH.LIBD.hippo.txt",quote = F,row.names = F,sep="\t")
            
sel=list.files(path = "E:/迅雷下载/20201222.LIBD.data.807ASH/rs",pattern="rs",full.name=T)
i=2	#因为第一个文件是空的
tmp=read.csv(sel[i],header = T)
tmp=dplyr::select(tmp,snpRsNum,feature_id,Type)
rt=tmp

for(i in 3:length(sel)){
  tmp=read.csv(sel[i],header = T)
  if(dim(tmp)[1]>0){
    tmp=dplyr::select(tmp,snpRsNum,feature_id,Type)
    rt=rbind(rt,tmp) 
  }
}

rsid=unique(rt$snpRsNum)
i=1
rt1=rt[rt$snpRsNum==rsid[i],]
feature.id=paste(unique(as.character(rt1$feature_id)),collapse = ":")
types=paste(unique(as.character(rt1$Type)),collapse = ":")
rt1=data.frame(avsnp150=rt1[1,1],feature.id=feature.id,types=types)
rts=rt1
for(i in 2:length(rsid)){
  rt1=rt[rt$snpRsNum==rsid[i],]
  feature.id=paste(unique(as.character(rt1$feature_id)),collapse = ":")
  types=paste(unique(as.character(rt1$Type)),collapse = ":")
  rt1=data.frame(avsnp150=rt1[1,1],feature.id=feature.id,types=types)
  rts=rbind(rts,rt1)
}
names(rts)=c("avsnp150","LIBD.eQTL.feature.id","LIBD.eQTL.type")
file1=merge(file,rts,by="avsnp150",all.x = T)
write.csv(file1,"at.least.two.ASH.in.DC.BF.10.add.eQTL.GWAS.add.DEG.LIBDGWAS.csv")
